% This function computes the conditional standard deviations of y at x_t using simulations
function out = conStdY_sim(xFull,model,pruningOn)

rng(1,'twister')
numSim    = 250;
shocks    = randn(model.ne,numSim);
shocks    = (shocks-mean(shocks,2))./repmat(std(shocks,0,2),1,numSim);


T         = size(xFull,2);
selectxf  = ones(model.nx,1);
selectnx1 = zeros(model.nx,1);
selectnx1(1:model.nx1,1) = 1;
if model.orderApp == 3
    xSim      = xFull([selectxf;selectnx1;selectnx1]==1,:);
elseif model.orderApp == 2
    xSim      = xFull([selectxf;selectnx1]==1,:);
elseif model.orderApp == 1
    xSim      = xFull([selectxf]==1,:);
end
    
if pruningOn == 1
    [~,y0]= modelFit_pruned(xSim(:,1),model);
else
    [~,y0]= modelFit(xSim(:,1),model);
end
ny      = length(y0);
stdy    = zeros(ny,T);
y_cu    = zeros(ny,T);
y_cup   = zeros(ny,numSim);
for t=1:T
    x_cu = xSim(:,t);
    if pruningOn == 1
        [~,y_cu(:,t)] = modelFit_pruned(x_cu,model);
        x_cup = stateTransition_pruned(x_cu,model);
    else
        [~,y_cu(:,t)] = modelFit(x_cu,model);
        x_cup = stateTransition(x_cu,model);
    end
    for s=1:numSim
        x_cup_sim = x_cup;
        x_cup_sim(1:model.nx,1) = x_cup(1:model.nx,1) + model.eta*shocks(:,s);
        if pruningOn == 1
            [~,y_cup(:,s)]   = modelFit_pruned(x_cup_sim,model);
        else
            [~,y_cup(:,s)]   = modelFit(x_cup_sim,model);
        end
    end
    stdy(:,t) = std(y_cup,0,2);
end

%% Output
out.ySim            = y_cu;
out.conStdY         = stdy;
labelYields = {};
for i=1:length(model.maturities)
    labelYields{i} = ['$r^{(',num2str(model.maturities(1,i)),')}_t$'];
end
out.labely  = [model.labely,labelYields];

end